In [1]:
import importlib
import theano.tensor as T
import sys, os
sys.path.append("/home/bl3/PycharmProjects/GeMpy/")
sys.path.append("/home/bl3/PycharmProjects/pygeomod/pygeomod")
import GeoMig
#import geogrid
#importlib.reload(GeoMig)
importlib.reload(GeoMig)
import numpy as np
os.environ['CUDA_LAUNCH_BLOCKING'] = '1'
np.set_printoptions(precision = 15, linewidth= 300, suppress = True)
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
In [ ]:
In [2]:
test = GeoMig.GeoMigSim_pro2(c_o = np.float32(-0.1),range = 17)
test.create_regular_grid_3D(0,10,0,10,0,10,20,20,20)
test.theano_set_3D_nugget_degree0()
In [3]:
layer_1 = np.array([[1,5,7], [9,5,7]], dtype = "float32")
layer_2 = np.array([[2,5,1],[7,5,1]], dtype = "float32")
dip_pos_1 = np.array([2,5,6], dtype = "float32")
dip_pos_2 = np.array([6.,4,6], dtype = "float32")
dip_pos_3 = np.array([8,4,5], dtype = "float32")
dip_angle_1 = float(0)
dip_angle_2 = float(45)
layers = np.asarray([layer_1,layer_2])
dips = np.asarray([dip_pos_1])#, dip_pos_3])
dips_angles = np.asarray([dip_angle_1], dtype="float32")
azimuths = np.asarray([0], dtype="float32")
polarity = np.asarray([1], dtype="float32")
#print (dips_angles)
rest = np.vstack((i[1:] for i in layers))
ref = np.vstack((np.tile(i[0],(np.shape(i)[0]-1,1)) for i in layers))
dips_angles.dtype
rest = rest.astype("float32")
ref = ref.astype("float32")
dips = dips.astype("float32")
dips_angles = dips_angles.astype("float32")
type(dips_angles)
Out[3]:
In [4]:
rest, ref
Out[4]:
In [5]:
rest, ref
Out[5]:
In [6]:
G_x = np.sin(np.deg2rad(dips_angles)) * np.sin(np.deg2rad(azimuths)) * polarity
G_y = np.sin(np.deg2rad(dips_angles)) * np.cos(np.deg2rad(azimuths)) * polarity
G_z = np.cos(np.deg2rad(dips_angles)) * polarity
G_x, G_y, G_z
Out[6]:
In [7]:
test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0]
In [ ]:
_,h1 = np.argmin((abs(test.grid - ref[0])).sum(1)), test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0][np.argmin((abs(test.grid - ref[0])).sum(1))]
_, h2 =np.argmin((abs(test.grid - ref[1])).sum(1)), test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0][np.argmin((abs(test.grid - ref[1])).sum(1))]
In [ ]:
print(test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0][np.argmin((abs(test.grid - ref[0])).sum(1))])
for i in range(rest.shape[0]):
print(test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0][np.argmin((abs(test.grid - rest[i])).sum(1))])
In [39]:
test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0][np.argmin((abs(test.grid - rest[0])).sum(1))]
Out[39]:
In [40]:
rest
Out[40]:
In [41]:
sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0].reshape(200,200,200, order = "C")[:,:,::-1].transpose()
#sol = np.swapaxes(sol,0,1)
In [42]:
G_x, G_y, G_z = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[-3:]
G_x, G_y, G_z
Out[42]:
In [43]:
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import matplotlib
import matplotlib.cm as cmx
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
h = np.array([h1,h2])
cm = plt.get_cmap("jet")
cNorm = matplotlib.colors.Normalize(vmin=h.min(), vmax=h.max())
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)
sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0].reshape(200,200,200,
order = "C")[:,:,:]
#sol = np.swapaxes(sol,0,1)
from skimage import measure
isolines = np.linspace(h1,h2,2)
#vertices = measure.marching_cubes(sol, isolines[0], spacing = (0.2,0.2,0.2),
# gradient_direction = "descent")[0]
for i in isolines[0:10]:
vertices = measure.marching_cubes(sol, i, spacing = (0.05,0.05,0.05),
gradient_direction = "ascent")[0]
ax.scatter(vertices[::40,0],vertices[::40,1],vertices[::40,2],color=scalarMap.to_rgba(i),
alpha = 0.2) #color=scalarMap.to_rgba(vertices[::10,2])
ax.scatter(layers[0][:,0],layers[0][:,1],layers[0][:,2], s = 50, c = "r" )
ax.scatter(layers[1][:,0],layers[1][:,1],layers[1][:,2], s = 50, c = "g" )
ax.quiver3D(dips[:,0],dips[:,1],dips[:,2], G_x,G_y,G_z, pivot = "tail", linewidths = 2)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.set_xlim(0,10)
ax.set_ylim(0,10)
ax.set_zlim(0,10)
#ax.scatter(simplices[:,0],simplices[:,1],simplices[:,2])
Out[43]:
In [13]:
test.c_o.set_value(-0.56)
test.c_o.get_value()
Out[13]:
In [81]:
test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[1]
Out[81]:
In [13]:
c_sol = np.array(([-7.2386541205560206435620784759521484375E-14],
[-1.5265566588595902430824935436248779296875E-14],
[-1.154631945610162802040576934814453125E-14],
[6.21724893790087662637233734130859375E-15],
[-5.9952043329758453182876110076904296875E-15],
[7.99360577730112709105014801025390625E-15],
[2.220446049250313080847263336181640625E-15],
[-3.641531520770513452589511871337890625E-14],
[8.0380146982861333526670932769775390625E-14],
[0.8816416857576581111999303175252862274646759033203125],
[9.355249580684368737593104015104472637176513671875],
[-0.1793850547262900996248191631821100600063800811767578125],
[0.047149729032205163481439313954979297704994678497314453125],
[-8.994519501910499315044944523833692073822021484375],
[ 0.4451793036427798000431721447966992855072021484375],
[-1.7549816402777651536126768405665643513202667236328125],
[0.0920938443689063301889063950511626899242401123046875],
[0.36837537747562587586713789278292097151279449462890625])).squeeze()
In [71]:
c_sol.squeeze()
Out[71]:
In [8]:
In [4]:
c_sol=np.array([ -0.07519608514102089913411219868066837079823017120361328125,
0,
3.33264951481644633446421721600927412509918212890625,
1.3778510792932487927231477442546747624874114990234375,
-2.295940519242440469582788864499889314174652099609375,
])
In [ ]:
In [5]:
import pymc as pm
a = pm.Uniform('a', lower=-1.1, upper=1.1, )
b = pm.Uniform('b', lower=-1.1, upper=1.1, )
c = pm.Uniform('c', lower=-1.1, upper=1.1, )
d = pm.Uniform('d', lower=-1.1, upper=1.1, )
e = pm.Uniform('e', lower=-1.1, upper=1.1, )
f = pm.Uniform('f', lower=-1.1, upper=1.1, )
@pm.deterministic
def this(value = 0, a = a ,b = b,c = c,d = d,e= e,f =f, c_sol = c_sol):
sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,
a,b,-3*b,d,e,f)[1]
#error = abs(sol-c_sol)
#print (error)
return sol
like= pm.Normal("likelihood", this, 1./np.square(0.0000000000000000000000000000000000000000000000001),
value = c_sol, observed = True, size = len(c_sol)
)
model = pm.Model([a,b,c,d,e,f, like])
In [6]:
M = pm.MAP(model)
M.fit()
In [8]:
a.value, b.value, c.value,d.value, e.value, f.value, this.value, c_sol, test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[1]
Out[8]:
In [72]:
a.value, b.value, c.value,d.value, e.value, f.value, this.value, c_sol, test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[1]
Out[72]:
In [40]:
a.value, b.value, c.value,d.value, e.value, f.value, this.value, c_sol, test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[1]
Out[40]:
In [457]:
a.value, b.value, c.value,d.value, e.value, f.value, this.value, c_sol, test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,a,b,-3*b,d,1,1)[1]
Out[457]:
In [459]:
-3*b.value
Out[459]:
In [406]:
a.value, b.value, c.value,d.value, e.value, f.value, this.value, c_sol, test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[1]
Out[406]:
In [372]:
a.value, b.value, c.value,d.value, e.value, f.value, this.value, c_sol, test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[1]
Out[372]:
In [368]:
a.value, b.value, c.value,d.value, e.value, f.value, this.value, c_sol, test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[1]
Out[368]:
In [363]:
a.value, b.value, c.value,d.value, e.value, f.value, this.value, c_sol, test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[1]
Out[363]:
In [ ]:
In [ ]:
In [192]:
test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,0,0,-0.33,0,1,1)[1]
Out[192]:
In [ ]:
In [ ]:
In [450]:
a.value, b.value, c.value,d.value,e.value,f.value, this.value, c_sol, test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,a,b,1,1,1,1)[1]
Out[450]:
In [74]:
a.value, b.value, c.value,d.value,e.value,f.value, this.value, c_sol, test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[1]
Out[74]:
In [19]:
importlib.reload(GeoMig)
test = GeoMig.GeoMigSim_pro2(c_o = np.float32(-0.1),range = 17)
test.create_regular_grid_3D(0,10,0,10,0,10,20,20,20)
test.theano_set_3D_nugget_degree0()
In [19]:
a.value, b.value, c.value,d.value,e.value,f.value
Out[19]:
In [26]:
import matplotlib.pyplot as plt
%matplotlib inline
G_x, G_y, G_z = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[-3:]
sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,a,b,-3*b,d,1,-1)[0].reshape(20,20,20)
def plot_this_crap(direction):
fig = plt.figure()
ax = fig.add_subplot(111)
if direction == "x":
plt.arrow(dip_pos_1[1],dip_pos_1[2], dip_pos_1_v[1]-dip_pos_1[1],
dip_pos_1_v[2]-dip_pos_1[2], head_width = 0.2)
plt.arrow(dip_pos_2[1],dip_pos_2[2],dip_pos_2_v[1]-dip_pos_2[1],
dip_pos_2_v[2]-dip_pos_2[2], head_width = 0.2)
plt.plot(layer_1[:,1],layer_1[:,2], "o")
plt.plot(layer_2[:,1],layer_2[:,2], "o")
plt.plot(layer_1[:,1],layer_1[:,2], )
plt.plot(layer_2[:,1],layer_2[:,2], )
plt.contour( sol[25,:,:] ,30,extent = (0,10,0,10) )
if direction == "y":
plt.quiver(dips[:,0],dips[:,2], G_x,G_z, pivot = "tail")
plt.plot(layer_1[:,0],layer_1[:,2], "o")
plt.plot(layer_2[:,0],layer_2[:,2], "o")
plt.plot(layer_1[:,0],layer_1[:,2], )
plt.plot(layer_2[:,0],layer_2[:,2], )
plt.contour( sol[:,10,:].T ,30,extent = (0,10,0,10) )
if direction == "z":
plt.arrow(dip_pos_1[0],dip_pos_1[1], dip_pos_1_v[0]-dip_pos_1[0],
dip_pos_1_v[1]-dip_pos_1[1], head_width = 0.2)
plt.arrow(dip_pos_2[0],dip_pos_2[1],dip_pos_2_v[0]-dip_pos_2[0],
dip_pos_2_v[1]-dip_pos_2[1], head_width = 0.2)
plt.plot(layer_1[:,0],layer_1[:,1], "o")
plt.plot(layer_2[:,0],layer_2[:,1], "o")
plt.plot(layer_1[:,0],layer_1[:,1], )
plt.plot(layer_2[:,0],layer_2[:,1], )
plt.contour( sol[:,:,25] ,30,extent = (0,10,0,10) )
#plt.colorbar()
#plt.xlim(0,10)
#plt.ylim(0,10)
plt.colorbar()
plt.title("GeoBulleter v 0.1")
In [40]:
sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,a,b,-3*b,d,+1,+1)[0].reshape(20,20,20)
In [41]:
plot_this_crap("y")
In [ ]:
In [436]:
a.value, b.value
Out[436]:
In [293]:
test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[1]
Out[293]:
In [440]:
c_sol
Out[440]:
In [ ]:
In [ ]:
In [ ]:
In [50]:
h,j,k =sol[5,10,35], sol[25,5,5], sol[30,15,-25]
layer_1 = np.array([[1,5,7],[5,5,7],[6,5,7], [9,5,7]], dtype = "float32")
layer_2 = np.array([[1,5,1],[5,5,1],[9,5,1]], dtype = "float32")
print(sol[5,25,35], sol[25,25,35], sol[30,25,35], sol[45,25,35])
print(sol[5,25,5], sol[25,25,5], sol[45,25,5])
In [51]:
list(layer_1[0]*5)
Out[51]:
In [52]:
interfaces_aux = test.geoMigueller(dips,dips_angles,azimuths,polarity,
rest, ref)[0]
h = sol[10,20,30]# interfaces_aux[np.argmin(abs((test.grid - ref[0]).sum(1)))]
k = sol[30,15,25]# interfaces_aux[np.argmin(abs((test.grid - dips[0]).sum(1)))]
j = sol[45,25,5]#interfaces_aux[np.argmin(abs((test.grid - dips[-1]).sum(1)))]
h,k,j
Out[52]:
In [53]:
dips[-1], ref[0]
Out[53]:
In [54]:
sol[30,15,25], sol[30,15,25]
Out[54]:
In [59]:
sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0].reshape(50,50,50, order = "C")
sol = np.swapaxes(sol,0,1)
plt.contour(sol[:,25,:].transpose())
Out[59]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [121]:
"""Export model to VTK
Export the geology blocks to VTK for visualisation of the entire 3-D model in an
external VTK viewer, e.g. Paraview.
..Note:: Requires pyevtk, available for free on: https://github.com/firedrakeproject/firedrake/tree/master/python/evtk
**Optional keywords**:
- *vtk_filename* = string : filename of VTK file (default: output_name)
- *data* = np.array : data array to export to VKT (default: entire block model)
"""
vtk_filename = "noddyFunct2"
extent_x = 10
extent_y = 10
extent_z = 10
delx = 0.2
dely = 0.2
delz = 0.2
from pyevtk.hl import gridToVTK
# Coordinates
x = np.arange(0, extent_x + 0.1*delx, delx, dtype='float64')
y = np.arange(0, extent_y + 0.1*dely, dely, dtype='float64')
z = np.arange(0, extent_z + 0.1*delz, delz, dtype='float64')
# self.block = np.swapaxes(self.block, 0, 2)
gridToVTK(vtk_filename, x, y, z, cellData = {"geology" : sol})
In [201]:
len(x)
Out[201]:
In [202]:
surf_eq.min()
Out[202]:
In [203]:
np.min(z)
Out[203]:
In [204]:
layers[0][:,0]
Out[204]:
In [ ]:
In [ ]:
In [ ]:
In [168]:
G_x = np.sin(np.deg2rad(dips_angles)) * np.sin(np.deg2rad(azimuths)) * polarity
G_y = np.sin(np.deg2rad(dips_angles)) * np.cos(np.deg2rad(azimuths)) * polarity
G_z = np.cos(np.deg2rad(dips_angles)) * polarity
In [183]:
a
Out[183]:
In [ ]:
data = [trace1, trace2]
layout = go.Layout(
xaxis=dict(
range=[2, 5]
),
yaxis=dict(
range=[2, 5]
)
)
fig = go.Figure(data=data, layout=layout)
In [ ]:
In [53]:
In [ ]:
import lxml
lxml??
In [11]:
# Random Box
#layers = [np.random.uniform(0,10,(10,2)) for i in range(100)]
#dips = np.random.uniform(0,10, (60,2))
#dips_angles = np.random.normal(90,10, 60)
#rest = (np.vstack((i[1:] for i in layers)))
#ref = np.vstack((np.tile(i[0],(np.shape(i)[0]-1,1)) for i in layers))
#rest;
In [ ]:
In [63]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y, Z = axes3d.get_test_data(0.05)
cset = ax.contour(X, Y, Z, cmap=cm.coolwarm)
ax.clabel(cset, fontsize=9, inline=1)
print(X)
plt.show()
In [11]:
import matplotlib.pyplot as plt
% matplotlib inline
plt.contour( sol.reshape(100,100) ,30,extent = (0,10,0,10) )
Out[11]:
In [131]:
import matplotlib.pyplot as plt
% matplotlib inline
dip_pos_1_v = np.array([np.cos(np.deg2rad(dip_angle_1))*1,
np.sin(np.deg2rad(dip_angle_1))]) + dip_pos_1
dip_pos_2_v = np.array([np.cos(np.deg2rad(dip_angle_2))*1,
np.sin(np.deg2rad(dip_angle_2))]) + dip_pos_2
plt.arrow(dip_pos_1[0],dip_pos_1[1], dip_pos_1_v[0]-dip_pos_1[0],
dip_pos_1_v[1]-dip_pos_1[1], head_width = 0.2)
plt.arrow(dip_pos_2[0],dip_pos_2[1],dip_pos_2_v[0]-dip_pos_2[0],
dip_pos_2_v[1]-dip_pos_2[1], head_width = 0.2)
plt.plot(layer_1[:,0],layer_1[:,1], "o")
plt.plot(layer_2[:,0],layer_2[:,1], "o")
plt.plot(layer_1[:,0],layer_1[:,1], )
plt.plot(layer_2[:,0],layer_2[:,1], )
plt.contour( sol.reshape(100,100) ,30,extent = (0,10,0,10) )
#plt.colorbar()
#plt.xlim(0,10)
#plt.ylim(0,10)
plt.title("GeoBulleter v 0.1")
print (dip_pos_1_v, dip_pos_2_v, layer_1)
In [17]:
%%timeit
sol = test.geoMigueller(dips,dips_angles,rest, ref)[0]
In [18]:
test.geoMigueller.profile.summary()
In [ ]:
In [23]:
sys.path.append("/home/bl3/anaconda3/lib/python3.5/site-packages/PyEVTK-1.0.0-py3.5.egg_FILES/pyevtk")
nx = 50
ny = 50
nz = 50
xmin = 1
ymin = 1
zmin = 1
grid = sol
var_name = "Geology"
#from evtk.hl import gridToVTK
import pyevtk
from pyevtk.hl import gridToVTK
# define coordinates
x = np.zeros(nx + 1)
y = np.zeros(ny + 1)
z = np.zeros(nz + 1)
x[1:] = np.cumsum(delx)
y[1:] = np.cumsum(dely)
z[1:] = np.cumsum(delz)
# plot in coordinates
x += xmin
y += ymin
z += zmin
print (len(x), x)
gridToVTK("GeoMigueller", x, y, z,
cellData = {var_name: grid})
Out[23]:
In [32]:
%%timeit
sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref);
In [33]:
test.geoMigueller.profile.summary()
In [30]:
importlib.reload(GeoMig)
test = GeoMig.GeoMigSim_pro2()
In [42]:
Out[42]:
In [17]:
In [1]:
from theano import function, config, shared, sandbox
import theano.tensor as T
import numpy
import time
vlen = 10 * 30 * 768 # 10 x #cores x # threads per core
iters = 1000
rng = numpy.random.RandomState(22)
x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
f = function([], T.exp(x))
print(f.maker.fgraph.toposort())
t0 = time.time()
for i in range(iters):
r = f()
t1 = time.time()
print("Looping %d times took %f seconds" % (iters, t1 - t0))
print("Result is %s" % (r,))
if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
print('Used the cpu')
else:
print('Used the gpu')
In [1]:
from theano import function, config, shared, sandbox
import theano.tensor as T
import numpy
import time
vlen = 10 * 30 * 768 # 10 x #cores x # threads per core
iters = 1000
rng = numpy.random.RandomState(22)
x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
f = function([], T.exp(x))
print(f.maker.fgraph.toposort())
t0 = time.time()
for i in range(iters):
r = f()
t1 = time.time()
print("Looping %d times took %f seconds" % (iters, t1 - t0))
print("Result is %s" % (r,))
if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
print('Used the cpu')
else:
print('Used the gpu')
In [18]:
from theano import function, config, shared, sandbox
import theano.tensor as T
import numpy
import time
vlen = 10 * 30 * 768 # 10 x #cores x # threads per core
iters = 1000
rng = numpy.random.RandomState(22)
x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
f = function([], T.exp(x))
print(f.maker.fgraph.toposort())
t0 = time.time()
for i in range(iters):
r = f()
t1 = time.time()
print("Looping %d times took %f seconds" % (iters, t1 - t0))
print("Result is %s" % (r,))
if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
print('Used the cpu')
else:
print('Used the gpu')
In [ ]:
In [759]:
np.set_printoptions(precision=2)
test.geoMigueller(dips,dips_angles,rest, ref)[1]
Out[759]:
In [751]:
T.fill_diagonal?
In [758]:
import matplotlib.pyplot as plt
% matplotlib inline
dip_pos_1_v = np.array([np.cos(np.deg2rad(dip_angle_1))*1,
np.sin(np.deg2rad(dip_angle_1))]) + dip_pos_1
dip_pos_2_v = np.array([np.cos(np.deg2rad(dip_angle_2))*1,
np.sin(np.deg2rad(dip_angle_2))]) + dip_pos_2
plt.arrow(dip_pos_1[0],dip_pos_1[1], dip_pos_1_v[0]-dip_pos_1[0],
dip_pos_1_v[1]-dip_pos_1[1], head_width = 0.2)
plt.arrow(dip_pos_2[0],dip_pos_2[1],dip_pos_2_v[0]-dip_pos_2[0],
dip_pos_2_v[1]-dip_pos_2[1], head_width = 0.2)
plt.plot(layer_1[:,0],layer_1[:,1], "o")
plt.plot(layer_2[:,0],layer_2[:,1], "o")
plt.plot(layer_1[:,0],layer_1[:,1], )
plt.plot(layer_2[:,0],layer_2[:,1], )
plt.contour( sol.reshape(50,50) ,30,extent = (0,10,0,10) )
#plt.colorbar()
#plt.xlim(0,10)
#plt.ylim(0,10)
plt.title("GeoBulleter v 0.1")
print (dip_pos_1_v, dip_pos_2_v, layer_1)
In [443]:
n = 10
#a = T.horizontal_stack(T.vertical_stack(T.ones(n),T.zeros(n)), T.vertical_stack(T.zeros(n), T.ones(n)))
a = T.zeros(n)
print (a.eval())
#U_G = T.horizontal_stack(([T.ones(n),T.zeros(n)],[T.zeros(n),T.ones(n)]))
In [6]:
T.stack?ö+aeg
In [13]:
x_min = 0
x_max = 10
y_min = 0
y_max = 10
z_min = 0
z_max = 10
nx = 2
ny = 2
nz = 2
g = np.meshgrid(
np.linspace(x_min, x_max, nx, dtype="float32"),
np.linspace(y_min, y_max, ny, dtype="float32"),
np.linspace(z_min, z_max, nz, dtype="float32"), indexing="ij"
)
np.vstack(map(np.ravel, g)).T.astype("float32")
Out[13]:
In [18]:
map(np.ravel, g)
Out[18]:
In [22]:
np.ravel(g, order = "F")
Out[22]:
In [23]:
g
Out[23]:
In [25]:
np.transpose?
In [117]:
from scipy.optimize import basinhopping
In [ ]:
c_sol, test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,1,1,1,1,1,1)[1]
In [127]:
def func2d(x):
return abs((test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref,x[0],x[1],x[2],x[3],1,1)[1] - c_sol)).sum()
In [128]:
minimizer_kwargs = {"method": "BFGS"}
x0 = [0.1, 0.1,0.1,0.1]
ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs,
niter=200)
In [123]:
ret
Out[123]:
In [126]:
ret
Out[126]:
In [129]:
ret
Out[129]:
In [ ]: